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ABSTRACT 


The present study aims to explore the possibility of employing the fast emerging Virtual 
Instrumentation (VI) technology for condition monitoring applications in rotatmg 
machinery. Time Series Analysis is known to be an effective mathematical tool, for 
analysis and forecasting of the behavior of dynamic systems. The objective of the present 
study is to translate known mathematical tools for time series analysis into Virtual 
Instrument by employing Graphical Programming Software, LAbVTEW The time series 
analysis procedure and the graphical programming features have been reviewed, briefly. 
Virtual Instruments have been developed to analyze standard sets of time-series data. The 
Vis developed, have been described in detail along with their various substructures The 
Vis presently incorporate Model Identification. Model Estimation and Forecastmg. 
Features like Diagnostic checking and Noise considerations could not be incorporated into 
the package during the present smdy. These can be built as the next stage of 
development, to make the package complete, for on-line usage on rotating machinery 
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CHAPTER 1 


INTRODUCTION 

The idea of using a mathematical method to describe the behaviour of a physical 
phenomenon is well estabhshed. It is sometimes possible to derive a model based on 
physical laws, which enables us to calculate the value of some time dependent quantity 
nearly exactly, at any instant of time. If exact calculation were possible, such a model 
would be entirely deterministic However, no phenomenon is totally deterministic and in a 
variety of problems, there are many unknown factors for which it is not possible to write a 
deterministic model that allows exact calculation of future behaviour of the phenomenon. 
In such cases, it may be possible to derive a model that can be used to calculate the 
probabihty of a future value lying between two specified values limits. Such a model is 
called a probabihty model or a stochastic model. 

The aim of the present smdy is to explore the potential offered by the fast emerging 
Virtual Instrumentation Technology for time series analysis Virtual Instrumentation offers 
the possibhity of a more eflacient and flexible mstrumentation. data acquisition and 
analysis, which can be very useful for condition monitoring apphcations, particularly, 
vibration based health diagnosis of rotor-bearing systems. 

Condition Monitoring and Health Diagnosis of rotating machinery in aero apphcations is 
an important area of technology development. Reference can be made to the texts by 
Childs (1993) and Rao (1992), for a discussion on the various aspects of rotating 
machinery modeling and health evaluation. The topics of mechanical fault monitormg and 
diagnosis have been comprehensively discussed by Mitcheh ( 1981) and CoUacott ( 1987). 

The progra mmin g flexibuity and the cost effectweness of the Virtual Instrumentation 
technology, offers a possibhity of more efficient and powerful vibration monitoring of 
rotor and forecasting of then dynamics. A traditional vibration analysis instrument is self 
contained, with signal input/output capabilities and fixed user interface features such as 
knobs, pushbuttons and other features. Specialized circuitry, including A/D converters, 
signal conditioning, microprocessors, memory and internal real-world signals, analyze 
them and present results to the user (e.g. oscilloscope). The user cannot change the 
instrument functionality Consider a typical rotordynamic application of forecasting the 
dynamics at time r+Ar on the basis of a vibration signature measured at time t The 
procedure involves transferring the signal measured on a traditional oscilloscope or signal 
analyzer, onto a computer in digitized form, writing a computer program, obtaining the 
required output in numerical form and finally using a graphics package to display the 
forecasted dynamics. On the other hand, graphical programming in conjunction with 
appropriate data acquisition cards, offers the possibhity of packaging the entire process 
above into a less expansive and more powerful virtual instrument simultaneously 
displaying the current and forecasted machine signature on-line and simultaneously on the 
computer screen. 
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The salient feature of a virtual instrument are 


• user defined 

• application oriented system with connectivity to networks, peripherals, and 
appHcations 

• software based 

• low cost, reusable 

• open, flexible functionahty based on familiar computer technology 

Figures 1.1(a) and 1. 1(b) depict the comparison between a traditional hardware instrument 
configuration and a Virtual Instrument configuration. 


Trajdifiona] Instnrn«ni 
COsciJloscop*) 



Rotor Di5c 


Fig. 1.1(a) Traditional Hardware Instrument configimation 



Rotor Di5C 


Fig. 1.1(b) Typical Virtual Instrument Configimation 
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The present study aims to bmld Virtual Instruments for typical Time Series Analysis 
applications, as in rotors. 

Time Series Analysis involves 

(i) Model Identification 

(ii) Parameter Estimation 

(iii) Diagnostic Checking 

(iv) Forecasting 

Fig. 1.2 outlines the procedure for Time Series Analysis. 



Fig. 1.2 Iterative Model Building Approach 


Chapter 2 gives a brief review of Time Series Analysis procedures. The Virtual 
Instrumentation developed during the course of this study is described in Chapter 3. The 
conclusions and scope for future work is outlined in Chapter 4. 
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CHAPTER 2 


TIME SERIES ANALYSIS: A REVIEW 


The use at time t of available observations from a series in time, to forecast its value at 
some future time, t+l, finds applications in a variety of engineering, business, economic, 
production planning. The time period, I, over which the forecast is required is called the 
'lead time'. For exanqjle. in a condition monitoring problem, the vibratory displacement 

reading r, at the present hour and the readings Zi. 2 , -t-i , m the previous hotirs 

might be used to forecast the vibratory displacement of the machine for lead times I ~ 
1,2,3 ..hours ahead. The accuracy of the forecasts is expressed by calculating the 
probability limits on either side of each forecast. These limit s are calculated for a 
convenient set of probabihties, generally between 50% to 90%. 

Tlie methods for obtaining the forecasts and esti m atmg their probabihty limits are 
described in this chapter. A stochastic process is generally usefully defined in terms of its 
mean, variance and autocorrelation function. These parameters are extensively used in the 
present study and for the sake of completeness, their definitions are provided below. 


Mean: The mean,// , of a series of N observations Zi, Z 2 , r\ recorded at times ti, 

U, ....iv is the expected value £'[Zf] ofrf 

and is estimated as 


1 ^ 

Autocovariance: Auto covariance at lag k is defibaed by 
Yk = = E[{z, - p){z^^k “ /^)] 

and is estimated as 

ZV'r --"X k = 0,1,2 ,K 


( 2 . 1 ) 


( 2 . 2 ) 


N ^ 

f=i 

Autocorrelation: Autocorrelation at lag k is given by 

* (2.3) 

and its estimated value is given by , 

Partial Autocorrelation function: It is the last coefBcient of an autoregressive process of 
order / and is estimated as 

for / = 1 
for l = 2,3,....L 


hi = 


1 


/-I 


/-I 




where 
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~ 1 ^ - 1 ^ 2 , , 1-1 


( 2 . 5 ) 


2.1 REGRESSIVE PROCESSES 
A linear model 


:: = +^2^2**" +a 

relating a dependent variable r to a set of independent variables x, pins an error term a, is 
often referred to as a regression model, and r is said to be regressed on x,. 

2.1.1 Autoregressive Model 

In this model, the variable is regressed on previous values of itself and the current value is 
expressed as a finite linear aggregate of previous values of the process and a shock a, . If 

the values of the process at equally spaced times tj-\J-l. are 

Zf, Zj_i, Zf_ 2 , and 3) are the deviations fi:om mean // (e.g. 

Zj = z, - p) then the autoregressive model is given by 

% = +^2^T-2'^ ( 2 - 6 ) 

The above AR model is written in concise form as 

( 2 . 7 ) 

where the AR operator <f){E) is 

= -(jipBP ( 2 . 8 ) 

B is the backward operator defined as 

B% = 3 ;., 

Equation (2.6) is also expressed as 
- = y/{B)a^ 
where 

= ( 2 . 11 ) 


( 2 . 9 ) 

( 2 . 10 ) 
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This is called an autoregressive (ARJ process of order p, in brief AR{p), containing 
p +2 unJknown parameters p\ p-, (j~^ which have to be estimated from data. 

The additional parameter cx^ is the variance of the white noise process a, . 

Autocorrelation function satisfies the difference equation 

= 0 

Le. 

Pk - 4'2Pk-2^---'^4' pPk-p fc > 0 • (2.12) 

2.1.2 Moving Average Model 

In this model the current value of the process is expressed as a finite linear aggregate of 
current and previous shocks a's 

=r=a,- (2-13) 

or concisely as 

P^eiB^a, (2.14) 

where diE) is the MA operator given by 

9{E) = -OqBP (2.15) 

The above is called a moving average process of order q , in brief MA( 9 ), containing 
q + 2 unknown parameters p', 8^, ,6^; <j~^. 

2.1.3 Mixed Autoregressive-Moving Average (ARMA) Model 

To achieve greater flexibility in fitting of actual time series, it is at times advantageous to 
include both AR and MA terms in the model. Thus 

-t = + <t>2%~2'^ ~ ^2^t-2~ (2-16) 

which can be expressed as 

<f>iB)^,=0{B)a^ (2.17) 
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This is called Autoregreesrve Moving average (ARMA) process of order {p,q). The 

model contains p+q + 2 irnknovi/’n parameters (9j, ,6^; al that are 

estimated jfrom data. 

2.1.4 General Autoregressive Integrated Moving Average (ARIIVIA) Model 

This model is generally employed to analyze processes which are non- stationary in time 
and is written as 

. (2 18) 

where 

V'p=(l-Byz, (2.19) 

The above is called an ARIMA process which is equivalent to stationary 

ARMAC/?,!?) model 

(/>(B)w, = d{B)a, (2.20) 

where 

( 2 . 21 ) 

i.e. iv, series is the d th difference of the original time series deviations 3) . 

2.1.5 Time Series Analysis Procedure 

The various steps in time series analysis can be listed as 

i) Postulation of general class of models 

ii) Tentative identification of the model and initial parameter estimation 

iii) Final parameter estimation 

iv) Use of model so estimated for forecasting. 

A summary of the properties of autoregressive, moving average, mixed ARMA processes 
is given ia Table 2. 1. 
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Table 2.1 


Summary of the properties of AR, iVIA, and Mixed ARMA processes 



Autoregressive 

Moving Average 

Mixed ARMA 

Model in terms 
of previous z ' s 

~^{B)z^ = 

e-^\B)z,=a, 

e-'iBmB)^ =a, 

Model in terms 
of previous a' s 

z=r\B)a, 

z = 6{B)a, 

z, = r\B)e(B)a, 

n weights 

finite series 

mfinite series 

mfinite series 

^ weights 

infinite series 

finite series 

infinite series 

Stationarity 

roots of (j){B) = 0 lie 

always stationary 

roots of (f){B) = OHe 

condition 

outside the unit circle 


outside the unit 
circle 

Invertibility 

always invertible 

roots of 9{B') = Olie 

roots of 6{B) = 0 He 

condition 


outside the unit circle 

outside the imit 
circle 

Autocorrelation 

infinite (damped 

finite 

infinite (damped 

function 

exponentials and / or 
damped sine waves) 


exponentials and / or 
damped sine waves) 
after first q-p lags 


tails off 

cuts off at lag q 

tails ofif 

Partial auto- 

finite 

infinite (dominated 

infinite (dominated 

correlation 


by damped 

by damped 

function 


exponentials and /or 
sine waves) 

exponentials and /or 
sine waves after first 


cuts off 

tails off 

q-p lags) 
tails off 


2.2 MODEL IDENTIFICATION 

Identification methods are rough procedures applied to a set of data to indicate the kind 
of representational model which can be further investigated. The aim of the identification 
is to obtain some idea of the values of p,d and q needed to identify the general linear 
ARI]VIA(jy,if,( 7 ) model to obtain initial guesses for parameters. 

A general ARIMA(/7,<i,^) model is expressed in equation (2. 17) as 

= e(B)a, 

The above equation can be written as 

with 
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( 2 . 22 ) 


0„=Ma-£4.,) 

)=1 

The following approach is used for identification 

• difference z, as many times as is needed to produce stationarity, hopefully reducing 
the process under study to mixed ARMA process 

(j)(5)w^ =0Q+e(5)a^ (2.23) 

where 

• identify the resulting ARMA process 

The above steps are carried out through use of autocorrelation functions(acf) and partial 
autocorrelation functions(pacf). These functions are also used at the estimation stage to 
provide starting values for iterative procedures employed at that stage. 

2.2.1 Identification of Process Stationarity 

The degree of differencing d, is identified by noting that, for a stationary mixed ARMA 
process of order {p,0,q), 

if{B)zt =QiB)at 

the acf satisfies the difference equation 
(t)(5)p^ =0 for A: 

Also, if 

KB)=nf^l (1-G,B) 


then solution of this difference equation for the is of the form 

p = A-^G^ + A 2 G 2 + +ApGp for k>q-p (2.24) 

If the process is stationary, the requirement that the zeroes of (j)(5) lie outside the unit 
circle (ref: Table 2.1), implies that the roots Gi,G 2 ,. lie inside the unit circle. 
Equation (2.24) shows that, in case of a stationary model in which none of the roots lie 
close to the boundary of unit circle, the acf will quickly die out for moderate and large k. 
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If a single real root, say G, , approaches unity, then 

G| = 1 - 6 for (5 « 1) and also positive for large k . 
Therefore from (2.24), autocorrelation at lag k is 


p^«.4l(l-^5) (2.25) 

It can also be noted from Eqn.(2.25) that the acf does not die out quickly but falls off 
slowly and very nearly linearly. Therefore, a tendency for the acf not to die out quickly is 
taken as an indication that a root close to unity may exist. The estimated autocorrelation 
function tends to follow the behaviour of the theoretical autocorrelation function. 
Therefore, failure of the estimated acf to die out rapidly might logically suggest that the 
underlying stochastic process is to be treated as as non-stationary in , but possibly as 
stationary in Vz^, or in some higher difference. It is assumed that the degree of 
differencing d, necessary to achieve stationarity, has been reached when the acf of 

Wf = V^z^ dies out fairly quickly. In practice, d is either 0, 1, or 2 and is usually 
sufficient to inspect the first 20 or so estimated autocorrelations of the original series and 
of its first and second difference. 

2.2.2 Stationary ARMA Process Order Identification 

For an autoregressive process of order p , AR(/7) , whereas the autocorrelation function 
of an AR process tails off, its partial autocorrelation function has a cut-off after lag p . 
Conversely, the acf of a moving average process of order q has a cut-off after lag q , 
while its pacf tails off. If both the acf and pacf tail off, a mixed process is suggested. 
Furthermore, the acf for a mixed process, containing a /ith order autoregressive 
component and q order moving average component, is a mixture of exponential and 
damped sine waves after the first q - p lags. Conversely, the pacf for a mixed process is 
dominated by a mixture of exponentials and damped sine waves after first p — q lags. 

A summary of the behaviour of autocorrelation and partial autocorrelation functions is 
given in table 2.2. 


Table 2.2 Summary Of Behaviours Of ACF And PACF For Various ARMA Models 

(a) 


Order 

( 14 , 0 ) 

(04,1) 

Behaviour of 

decays exponentially 

only pj nonzero 

Behaviour of <j) 

only (j) j j nonzero 

exponential dominates 
decay 

prelim estimate from 

<t>i=Pi 

p, ^-ei/o+ej*) 

Admissible region 

-l<^l <1 

-i<ei <1 
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0 ?) 


Order 

(2,d,0) 

(0,d,2) 

Behaviour of 

mixture of exponentials or 
damped snie wave 

only p-^ nonzero 

Behaviour of (j) 

only (j)Yi and ^22 

nonzero 

dominated by mixture of 
exponentials or damped sine 
wave 

preUm estimate firom 

P^{l-P2) 

1 .2 

1-Pl 

. .Pi- Pi . 

I .2 

_ -^1(1- ^ 2 ) 
l + fff+0^ 

n - ~^2 

^ 1 + + ^2 

Admissible region 

-1<^z 52 <1 
^2 +^z^l < 1 
^2 “ ^1 ^ ^ 

-1<^2 
^2 1 
^2 -^1 <1 


(C) 


Order 

(l,d,l) 

Behaviour of p^ 

decays exponentially firom first lag 

Behaviour of 

dominated by exponential decay firom first 
lag 

prelim estimate firom 

Pi~ 0 P'i~ P\^\ 

l + df ~2pi0i 

Admissible region 

-l<^^j<l -1<0^<1 


2.2.3 Initial Parameter Estimation Of The ARMA Process 

Initial parameter estimation of an AJRMA {p.q) process is based on the first p + q + l 
autocovariances c j {j = 0,1 ,... , p + q) of . It is carried out in the following 

three stages: 

• AR parameters p estimated fi-om the autocovariances 

by sohdng the p linear equations 


A(^ = X 


where 

^q+i 


i,j = \ 2,....,p 


(2.26) 
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Using the estimates (j)' s obtained in the previous step, the first q + l autocovariances 

I 

Cj (7 = 0,1, . . . . , (y) of the derived series 

I 

Wf -<t>p-^t-p 

I 

are calcixlated by treatmg the process = (l>{B)w^ as a moving average one i.e 


Wi 


I I 

0 {B)af . To start with, the auto covariance Cj of is expressed in terms of 


the autocovariances c of w, , as 


p p 

c , = ZT P>^ = 

^ /=0A=0 

p = 0 


(2.27) 


= c 


j 


The autocovariances ^q ^sed in an iterative calculation to 

compute initial estimates of the moving average parameters 6 ^-, 0 2,^2^ ’^<7 

2 

of the residual variance cr^ , through Newton-Raphson technique, as follows: 
Denoting the set of iterative values 

r'={TQ,v^., 

1 

0' 


.Vq) 


2 2 

where rn = a a ^ j - 


0 


7 = U,. 




(2.28) 


if is the estimate of r , through the Newton-Raphson algorithm, obtained at rth 
iteration, the new \ alues at the (; + l)th iteration are obtained by 

where 


= -(r')-Vz 


(2.29) 


f = ifQ,fl,-,fq), 

q-J ^ 

fi= 


(2.30) 


/=0 
and 


i i+j J 



^0 ^1 • • • '<7-2 V-1 ^<7 


^0 ^1 ^2 • • • • '^‘7 


^1 ^2 ■ • ■ ^q-\ ^*7 ^ 


0 Tq ri . . . . T 2 

T = 

'^2 '^3 ■■■ ^q ® ^ 

+ 

0 

0 

0 


1 

• 0 

• 0 

• 0 

• 0 

• 0 

0 

• 0 

. 

1 

[ 

i 

0 0 0 . . . . Tq 


(2.31) 
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Iterations are carried out with starting values Tq = = ^2 =....= = 0 till 


convergence (Le. 


/ 


becomes negligible. 7 = 0, 1 ,....,^). 


The estimate for the overall constant, 6^ , Jfrom Eqn.(2.22), is 


= w 


w 


1-Z^A forp>0 
V /=1 J 

for p = 0 


(2.32) 


Also, the estimate of white noise variance is: 


^2 _ .,2 


for q > 0 


foT q = 0 

i=i 


(2.33) 


2.3 MODEL ESTIMATION 

The estimation theory is based on the concept of likelihood function. If a sample of' 
.V observations is r . with which an iV-dimensional random variable can be associated and 
^ is a set of p+q + l parameters, (^, (9, /i) , the hkelihood function gives the 

conditional probability of q for the observed set z . 

2.3.1 Conditional Likelihood For Am ARIMA Process 

For an ARIMA model (p,d,q), the N = n+d original observations^, form a time series 
denoted by ,ro,Cj,C 2 From these observations, a series w, of 

n = N-d differences is generated where =V^Zf.. 

The general problem of fitting parameters (j) and ^ of an ARIMA model is, then, 
equivalent to filttiag to the the stationary invertible ARMA(p,^) model which is 
written as 

=^t - p^t-p+d^a,_^ +d2a,_2+....+0^a,_^ (2.34) 

where 

E[w,] = p 
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For d>Q, }a = 0 is often a good approximation, otherwise p = ^ w, In. 

/=! 

The w's cannot be substituted immediately in (2.34) because of the difficulty of starting 

up the difference equation (2.34). However, if p values w. of and q values a. of 
a's , prior to the commencement of the w series are given, then the values of 

,a„ conditional on this choice can be calculated in turn from equation (2.34). 

For any given choice of parameters <j),0 and of the starting values w*, o* a set of 
values a,(<|),0iw»,£,,w,) can be calculated successively. If the a's are normally 
distributed, the joint probability density function is 

p{a^,a 2 ,...,a„) oc a^exp{-(Ea/ /2 ct„^)} 

t=\ 

Taking log on both sides of the above proportionality, the log likelihood associated with 
the parameter values ((jj.OjCr^) , for a particular set of dataw, conditional on the choice 
of (w*,q*) is obtained as 

S'*((|),0) 

I*(^,Q,oJ = -n]na^ ( 2 - 35 ) 

where ‘S'*((t),0) = j}f) (2.36) 

t=i 

It can be noted firom Eqn.(2.35) that lower the ^.((l),^ , higher will be the log likelihood 
function. Also data are involved only through 5. , the sum of squares. In order to get the 

n 

maximum likelihood function, the sum of the square function ^ is minimized. The 

r=l 

minimized S*, called least square estimate, corresponds to the desired set of parameters. 
2.3.2 Nonlinear Estimation 

For most situations, the maximum likelihood estimates are closely approximated by the 
least square estimates, which make 

s.»,e)= i = E[a,]' 

/=-00 /=-00 

a minimum. 

n 

In practice, the infinite sum can be replaced by a finite sum . 

t=-Q 
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(2.41) 


p q 

;=I 7=1 

• For the given values of the parameters (p.,(j),03 the residual sum of squares is 

calculated from • 

t=-Q 

• The least squares estimates are calculated next. The values of the parameters which 
minimize the residual sum of squares are obtained by a constrained optimization 
method (Marquardt method). Denoting all the initial assumed set of parameters 

by 

P = (Pi>p 2 .--” = Pi) where k = p + q + \ (2.42) 

the derivative x,, is calculated from Eqn.(2.39) 



- Fig. 2.1 Backforecasting w series 
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Minimization is carried out in three iterative stages: 


(1) With a, , X, , , the ky.k matrix A = {A^j } is calculated as 

n 

t—Q 

( 2 ) The vector g with elements g,,g 2 > ’8k is calculated as 

n 

8 , = 

t=-Q 

(3) Scaling quantities (for future use) are computed as 

A ~ VA/ 


(2 43) 


(2.44) 

(2.45) 


Staged) 

(1) The modified (scaled and constrained) linearized equations 


A* h* =g* (2.46) 

are formed according to 

A* =1 + 71 (2.47) 

8:=8JD, 

(2) The equations are solved for /z* which is scaled back to give the parameter 
corrections 

hj = h] / Dj (2.48) 

(3) The new parameter values are calculated using 

P = p^+A (2.49) 


which is used in computing the next sum of squares iS'(P) . 

Stager3i 

(1) If iS(P)<S'(P^) and the parameter corrections h are all negligible then 
convergence is assumed. Otherwise, next iteration is started from Stage(l), with 
the new parameter values P and the reduced value of tx by a factor 10 . 

(2) If 5'(P) > jS'(P q) , the parameter tz is increased by a factor 10 and computation 
is started at Stage(2). 
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2.4 FORECASTING 


A general ARIMA process can be forecast in three different explicit forms: 

1 . in terms of previous values of the z' s and current and previous values of the ds by 
direct use of the difference equation. 

2. in terms of current and previous shocks only. 

3. in terms of a weighted sum of previous values z^_^ of the process and the current 
shock . 

Forecasting is done in the present study using the difference equation form (1 above) , of 
the model, which is 

=9i^,_i+ ^p^d^t-p-d -6 A-i- (2-50) 

The value Z,,;, />1 is forecasted at the time origin t . This forecast is said to be made at 
origin t for lead time / , as follows 

z,+; =cpiz,^/_i+ +<i? -e +a,^i (2.51) 

2.4,1 Minimum Mean Square Error Forecasts 

The minimum mean square error forecast at origin t , for lead time / is the conditional 
expectation of z,+; at time i (Box and Jenkins, 1976). 

Mathematically, by regarding z, (/) as the forecast function of I for a fixed origin t , 

Zfil) = 

(2.52) 

= 9 1 {Zi+i_p_(i'\ - 6 1 [«/+/_! ]---0 q 1 



where 



II 

o 

to 


7= 1, 2 

1 

I! 

r — 1 

Cl 

I 

11 

1 

II 

! 

7 = 0, 1,2,... 

11 

+ 

11 

O 

7 = 1, 2, 
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For moving average terms in the difference equation, the forecasting process is started off 
initially by setting unknown a' s equal to their imconditional expected values of zero. In 
general, if the moving average operator B{B) is of degree q, the forecast equations for 

2 ,( 1 ), f,(2), f,(^) depend directly on the a's but forecast at longer lead times do 

not. 

2.4.2 Forecast Errors 


The forecast error for lead time / and its variance are 
et(J) = +¥ M^t+1 

F(/) = var[e^ (/)] = (! + 11/ 1 +'^ 2 '^ 


(2.54) 


The weights, ly /_! iii equations (2.54) are calculated by equating 

coefficients of powers of B in 

(l-9j5-...-(p^^^5^+^)(l+M/j5+H/25^+...) = l-ei5-e25^-...-e^5^ 
to get 

ij/ y = 1 y = 0 (2.55) 

J 

= E^cp/Vy_/-0y y>l 


2.4.3 Probability Limits Of The Forecasts 


Since the square of the error between the forecasted value for lead / and the actual value 
is 

If the a' s are normal, the conditional probability distribution p{zt+i I > ) of a 

future value of the process will be Normal with mean z, (J) and standard deviation 


V{1) = 


/-I 




y=i 


a 


2 

a 


The upper and lower probability limits are, therefore, obtained by (Box & Jenkins, 1976) 
^,+,i±} = l,(I)±uJV(f) (2.56) 
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where 


w = 0.68 for 50 % probability 
= 1.65 for 90 % probability 
= 1.96 for 95 % probability 
= 2-58 for 99 % probability 

The forecast is updated with the arrival of every new data , by updating the origin to 
/ + 1 and by calculating the new forecast error. 

The Virtual Instruments developed on the basis of the Time Series Analysis procedure 
described here are discussed in Chapter 3. 
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CHAPTER 3 


VIRTUAL mSTRUMENTATION 

The original plan, in the present study, was to develop Virtual Instruments for Model 
Identification, Parameter Estimation, Diagnostic Checkmg and Forecasting of the 
dynamics of an available laboratory rotor rig. However, due to time constraints the 
attempt had to be restricted to off-line development of Virtual Instrumentation for general 
Time Series Analysis. The VTs have been tested and validated, on the basis of available 
literature and are in a state, where they can be implemented on the rotor-rig vibration 
signals, with Httle effort. The rotor-rig, plug-in data acquisition card and the software are 
briefly described in the following. 

3.1 THE ROTOR RIG 

The available laboratory rotor rig comprises of a shaft supported in bearings (with 
provision for rolling element or fluid film beatings), with discs and flexible couplings.lt can 
be employed to simulate a variety of real life rotor dynamic problems (e.g. unbalance, oil 
whirl, shaft; bow, cracks, bearing wear and tear etc.), refer Figure 3.1. Vibration signals 
can be picked up through accelerometers (or proxinuty pick ups). Signal conditioning is 
done through Charge Amplifiers. The conditioned and amplified signals fi-om Charge 
Amplifiers can be taken into the available PC486, through AD/DA conversion cards. 
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3.2 PLUG EV DATA ACQUISITION 


A Plug-in DAQ board, AT-MIO-16E-10, from National Instruments, Texas, has been 
acquired. Tliis board has various combinations of analog, digital, and timing inputs and 
outputs Tlie board can acquire data in three modes: 1) continuous acquisition of a single 
channel, 2) multichannel acquisition with continuous sc annin g, or 3) multichannel 
acquisition with interval scanning. There are sixteen input channels and eight output 
channels 

3.3 SOFTWARE 

A graphical software. Lab VIEW has been acquired from National Instruments, Texas, 
USA. Lab VIEW is a program development appHcation, much Klee various commercial C 
or BASIC development systems While other programming systems use text based 
languages to create the code. Lab VIEW uses a graphical programming language, G, to 
create programs in block diagrams. It uses icons and graphical symbols to describe 
programming action. LabVTEW contains functions and subroutines for various 
programming tasks, like in other languages. Lab VIEW programs are called Virtual 
Instruments (Vis) because their appearance and operation imitate actual instruments. They 
are, however, analogous to functions from conventional language programs. Vis have 
both an interactive user interface and source code equivalent and accept parameters from 
higher level Vis The following are the descriptions of the salient VI features. 

• The interactive user interfece of a VI is called Front Panel, because it simulates the 
panel of a physical mstrument. The front panel can contain graphs, push buttons, 
knobs and other control indicators. Data can be input using a mouse or a keyboard. 

• VT receives instmetions from a block diagram which is a graphical program written in 
G The block diagram is also the source code. 

• Vis are hierarchical and modular and can be used as top-level programs or as 
subprograms within other programs or subprograms. A VI within a VI is called a 
sub VI. 

The variety of functions available m Lab VIEW are described below. 

Arithmetic functions : All types of arithmetic functions like add, subtract, multiply, 
round to nearest, increment, decrement, scale by, power of 2, quotient and remainder 
round to + infinity, round to - infinity, complex conjugate, absolute value, negate, square 
toot, reciprocal, random munber etc. 

Logical functions : AND, OR, EXCLUSIVE OR, NOT, NOR, NAND etc. 

Trignometric and logarithmic functions : Sine, Cosine, tangent, inverse sine etc. 
Comparison function : Equal to, not equal, greater than, lesser than, max and 
min, greater than 0, not equal to zero etc. 

Conversion functions : To byte integer, to word integer, to long integer, to unsigned 
byte, complex to polar, cluster to array etc. 

String function : String length, concatenate, string, string subset, split string , to 
decimate etc. 
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Array Functions : array size used to find the length of an anay, index array to extract 
an element at the given index, replace array element to replace the existing element with 
the desired element, array subset to extract subset of given length of rows and columns, 
build arrays to build an array of elements etc. 

Lab VIEW has two stmctures to repeat execution of a subdiagram - the While Loop and 
the For Loop 

While Loop • Tliis is a post-iterative test loop structure that repeats a section of code 
until a condition is met (Fig 3.2(a)). It is comparable to a Do loop or a Repeat-Until loop 
in conventional programming languages. 



Fig. 3.2(a) Block diagram for While loop 


The For Loop : This loop (Fig. 3.2(b)) executes the diagram inside its border a 

predetermined number of times. It has two terminals : the count terminal ( an input 

terminal) and the iteration terminal. The count tenninal specifies the number of times to 

execute the loop. The iteration terminal contains the number of completed iterations. The 

For Loop is equivalent to 

For i = 0 toN - 1 

Execute Diagram Inside The Loop 



Fig. 3.2(b) Block diagram for loop 
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LabVIEW has two structures that control data flow. They are the Case 
structure( analogous to‘if-theu-else’ statement in a conventional text based language) and 
the Sequence stmcture 

Case structure . Tliis is a conditional branching control structure (Fig. 3.2(c)), which 
executes one and only one condition of its subdiagrams based on its input. It is the 
combination of the IF, THEN, ELSE, and CASE statements in control flow languages. 




if (number.ge 0) then 
s q_root=sq rt( n u m b e r) 
else 

sq_root=sqrt(abs(numl: 

endif 


Fig. 3.2(c) Block diagram for Case structure 


Sequence structure ; This is a program control structure that executes its subdiagrams in 
a specific, numeric order. It is commonly used to force nodes that are not data- dependent 
to execute in a desired order. The portion of the diagram to be executed first is placed in 
frame 0 of the Sequence structure, the diagram to be executed second is placed in frame 1, 
and so on Frame is subdiagram of a Sequence structure and only one frame is visible at a 
time. 

Some other relevant features of the software are - 

Strings: A string is a collection of ASCII characters. A numeric data can be passed as 
character strings and then these are converted to numbers. 

File I/O : These are fimctions for working with files. In addition to reading and writing 
data, the LabVIEW file I/O functions move and rename files and directories, create 
spreadsheet-type files of readable ASCII text, and write data hr binary form for speed and 
compactness. 

Code Interface Node (CIN); CIN is a node which links external code written in a 
conventional programming language to LabVIEW. A CIN appears on the diagram as a set 
of input and output terminals. 
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3.4 COMPUTER PROGRAMS 

Virtual instruments have been developed for six sets of time series data available in 
literature (Box and Jenkius, 1976). Hiis is necessary to validate and check the graphical 
programs before their on-hne application to the rotor-rig. Tlie sLx time series' taken for this 
purpose pertain to 

Series A. Chemical process concentration readings, every' two hours 

Series B IBM common stock closing prices, daily, 17th May '6 1- 2nd November '62 

Series C: Chemical process temperature readings, every minute. 

Series D: Chemical process viscosity readings, every hour 

Series E. Wolfer sunspot numbers, yearly 

Series F: Series of 70 consecutive yields from a batch chemical process 

These sets of data have been chosen for programming and validation for they represent a 
wide spectrum of/a7cnm stochastic behaviour. Series A is a typical ARMA( 1,0,1) process. 
Series B is an Integrated Moving Average, IMA(0,1,1) process. Series C is an Integrated 
Autoregressive, IAR(1,1,0) process. Series' D,E and F are pure AR(1,0,0), AR(2,0,0), 
AR(2,0,0) processes respectively. 

The overall picture of the Vis developed can be seen from Figures 3.3-3. 8. These figures 
show the front panels of the Virtual Instruments developed, for Series A-F, respectively. 
The panels simultaneously show the input data, the autocorrelations, the partial 
autocorrelations, the forecast and the forecast limits. Some of the estimated parameters 
can also be seen in the form of pellets. The forecast are shown for a probability of 50%. 

Figures 3.9(a) and 3.9(b) typically show the graphical program (called the Block 
Diagram), in a condensed form, for the front panels of Fig. 3.3 (Series A). T his main VI 
contains several sub Vis. Figure 3.10 describes the position of these sub Vis in the 
hierarchy of computational operations. Figure 3.11 gives a consohdated hst of these 
sub Vis. 

The various sub Vis which are used in the main VI are described below. 

del.vi (Fig. 3. 12): differences the original time series (input data) 

r = {zi,Z 2 ,. ■,Zj^}to get w = according to the 

relation = V '^z, 

c & r.vi (Fig.3. 13) calculates the autocovariance and autocorrelation coej05cients 

according to Eqns.(2.2) and (2.3). 

pacf.vi represents the graphical program to calculate partial 

(Fig.3. 14) autocorrelation according to Eqns.(2.4)-(2 5). 
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autoreg para.vi 
(Fig.3.15) 

calculates initial autoiegressive parameters according to 

Eqns.(2 26). 

prelim parameter 

estimationl.vi 

(Fig.3.16) 

Figs 3 16(a)-(f) together \vith Fig.3.15 represent the graphical 
program for preliminaiy estimation of autoregressive and moving 
average parameters at the tentative identification stage according 
to Eqns.(2.27)-(2.31). 

phi to phi’^.vi 
(Fig.3.17) 

Figs.3. 17(a)-(b) represent the program to convert <j) to (p 
according to (p{B) = - B)^ 

forecastA.vi 
(Fig 3 18) 

calculates forecasts of time series at given origin t for lead / 
according to Eqns. (2.51)-(2 53) (used as a sub VI of Fig 3.19) 

back fore.vi 
(Fig 3.19) 

calculates backforecasts of w series before starting the recui'sive 
calculation of 

a(t) & sum of sq.vi 
(Fig 3.20) 

recursively calculates a, and its sum of squares S according to 
Eqns.(2.41) and (2.36) at the model estimation stage. 

a(t) & residual sum 
of squares.vi 
(Fig.3.21). 

This is a program to recursively calculate a, and its sum of 
squares S , according to Eqns.(2.41) and (2.36) at the model 
estimation stage. This VI uses Fig. 3. 19 and Fig. 3. 20 as sub Vis. 

final estimation 

icon.vi 

(Fig.3.22) 

Figs.3. 22(a)-(g) represent the graphical program in LabVIEW 
incorporating the model estimation algorithm described in Eqns. 
(2.38), (2.41)- (2.49). This program uses the sub Vis of Fig.3.21 

forecast.vi 
(Fig 3 23) 

Figs. 3.23(a)-(i) forecasts the final model according to the 
algorithm described in Eqns. (2.50)-(2.56). 
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3.5 RESULTS 


Tlie results obtained for the Series A, after implementation of the Vis described above are 
discussed here 

Model Identification of series A Tlie behaviour of autocorrelation and partial 
autocorrelation functions are studied, as discussed in Chapter 2, to identify a general 
ARIMA (j),d,q) model. Tlie autocorrelation and partial autocorrelation coefficients 
obtained fiom the Vis of Figs.3 13 and 3. 14 are given in Tables 3. 1 - 3 2 and Figs. 3.24 - 
3.25 ■ 


Table 3.1 Estimated Autocorrelations for series A 


Lag/C 

^=0 

d=l 

d=2 

1 

0.57 

-0.41 

-0.65 

2 

0 50 

0.02 

0.18 

3 

0.40 

-0.07 

-0.04 

4 

0.36 

-0.01 

0.04 

5 

0.33 

-0.07 

-0.04 

6 

0.35 

-0.02 

-0.04 

7 

0 39 

0.15 

0.13 

8 

0.32 

-0.07 

-0.11 

9 

0.30 

0.04 

0.04 

10 

0.26 

0.02 

0.02 

11 

0.19 

-0.05 

-0.02 

12 

0 16 

-0.06 

-0.02 

13 

0.20 

1 

O 

o 

-0.04 

14 

0.24 

0.16 

0.18 

15 

0.14 

-0.17 

-0.19 

16 

0 18 

0.03 

0.08 

17 

0 20 

0.01 

-0.03 

18 

0.20 

0.08 

0.09 

19 

0 14 

-0.12 

-0.17 

20 

0.18 

0.15 

0.20 
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Fig. 3.25(c) Partial autocorrelation for d = 2 



It is observed that the autocorrelation for cf = 1 are small after the first lag, suggesting that 
tlfis series might be described by an IMA(0,1,1) process. However, from the 
autocorrelation function for J = 0, it is seen that after lag 1, autocorrelation decreases 
fairly regularly. Tlierefore, an alternative possibility is that the series is a mixed ARMA of 
order (1,0,1). Tlie partial autocorrelation function for J = 0 tends to support tlfis 
possibility. Tlius, for further analysis of series A, these two possible models EMA (0,1,1) 
and ARMA( 1,0,1) are considered. 

Initial Parameter Estimation : With tentative models IMA (0,1,1) and ARMA( 1,0,1) 
identified, the initial estimates of 6 and (f) are calculated from Eqns.(2.25)-(2.34). The 
estimates are given in Table 3.3. 

Table 3.3 Initial Estimates For Series A 


Model 


e 


^0 

MHMM 

0.87 

0.48 


2.25 


- 

0.53 

0.107 

0 00 


Final Estimation : With tlie tentative model identification and initial estimates of 
parameters, more efiicient estimates are calculated for the series, by using the initial 
models and parameters as initial guesses, e.g. for IMA(0,1,1), the initial guess is 
where (9 = 0.53 and // = 0. - \ 
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Nonlinear least square estimates are, now calculated using constrained optimization 
(Marquardt Method) by Eqns.(2.38), (2.4])-(2 49) and the graphical program of 
Figs 3.21(a)-(g). Tlie final estimates are given in Table 3.4 


Table 3.4 Final Estimates For Series A 


Model 


e 


^0 

(1,0,1) 

0.91 

0.56 

0 099 

1.62 

(0,1,1) 

“ 

0.71 

0.102 

0.00 


Forecasting: 


Forecasts at origin 144 with lead /= 1,2,. ..,53 are calculated, for ARJV[A( 1,0,1) possibility, 
according to Eqn.(2 52) The values of upper and lower limit s for 50% probability are 
calculated from Eqns.(2.53)-(2.56). The actual data, along with the corresponding 
forecast and the upper and lower limits are given in Table 3.5. 

Table 3.5 Forecasting At Origin 144 For 50% Probability Limit 


Actual 

Data 

Forecast 

Upper 

Limit 

Lower 

Limit 

16.70 

16.73 

16.95 

16.52 

16.90 

16.77 

16.99 

16.54 

17.40 

16.79 

17.03 

16.56 

17.10 

16.82 

17.06 

16.57 

17.00 

16.84 

17.09 

16.59 

16.80 

16.86 

17.12 

16.61 


16.88 

17.14 

16.62 

17.20 

16.90 

17.16 


17.40 

16.91 

17.18 

16.65 

17.20 


17.19 

16.66 

16.90 


17.21 

16.67 

16.80 


17.22 

16.68 

17.00 


17.23 

16.69 

17.40 


17.24 

16.70 

17.20 


17.25 

16.71 



17.26 

16.71 



17.27 

16.72 



17.27 

16.73 



17.28 

16.73 



17.29 

16.74 



17.29 

• 16.74 

16.90 

17.02 

17.29 

16.75 
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16.90 

17.02 

17.30 

16.75 

17.00 

17.03 

17.30 

16.75 

16.70 

17.03 

17.31 

16.76 

16.90 

17 03 

17.31 

16.76 

17.30 

17 04 

17.31 

16.76 

17.80 

17.04 

17.31 

16.76 

17.80 

17 04 

17.32 

16.76 

17.60 

17.04 

17.32 

16.77 

17.50 

17 04 

17.32 

16.77 

17.00 

HlBSHi 

17.32 

16.77 

16.90 

IHDSEEHi 

17.32 

16.77 

17.10 


17.32 

16.77 

17.20 

17.05 

17.32 

16.77 

17.40 

17.05 

17.33 

16.77 


17.05 

17.33 

16.78 

HHEE9H 

17.05 

17.33 

16.78 



17.33 

16.78 

17.00 

17 05 

17.33 

16.78 

17.00 

17.05 

17.33 

16.78 

17.20 

17 05 

17.33 

16.78 

17.30 

17 06" 

17.33 

16.78 

17.40 

17.06 

17.33 

16.78 

17.40 

17.06 

17.33 

16.78 

17.00 

17.06 

17.33 

16.78 

18.00 

17 06 

17.33 

16.78 

18.20 

17.06 

17.33 

16.78 

17.60 

17.06 

17.33 

16.78 

17.80 

17.06 

17.33 

16.78 

17.70 

17 06 

17.33 

16.78 

17.20 

17.06 

17.33 

16.78 

17.40 

17.06 

17.33 

16.78 


Similarly the remaining series’ were analysed. Summary of models identified for series A- 
F in terms of various initial parameter estimates are given in Table 3.6. 

Final estimates for all 6 series’ are given in Table 3.7. 
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Table 3.6 Initial Estimates For Series A-F 


Series 

Model 


e 

^2 

^0 

A 

(1,0,1) 

(zJ, = 0 87 

(9, =0 48 

0.098 

2.25 

B 

(0,1,1) 

- 

(9i =-0 09 

52 151 

-0.28 

C 

(1,1,0) 

=081 

- 

0.019 

-0.01 

D 

(1,0,0) 

(j)^ = 0.86 

- 

0.093 

1.27 

E 

(2,0,0) 

<zSi = 132 
?i,=-063 

- 

289.214 

14 86 

F 

(2,0,0) 

= -032 

^>2=018 


114.719 

i 

58.29 


Table 3,7 Final Estimates For Series A-F 


Series 

Model 


e 


^0 

A 

(1,0,1) 

o 

II 

= 056 

0.099 

1.62 

B 

(0,1,1) 

- 

(9i = -0.09 

52.438 

-0.28 

C 

(1,1,0) 

= 0.82 

- 

0.018 

-0.01 

D 

(1,0,0) 

=0.87 

- 

0.091 1 

1.17 

E ^ 

(2,0,0) 

(j)^ = \Al 
<j)^ =-073 


234.867 

14.72 

F 

(2,0,0) 

<!}j = -035 
«52 = 019 


117.751 

59.05 


3.6 REMARKS 

The numerical results obtained above, for the six sets of time series data are compared 
with those given by Box and Jenldns (1976). A good agreement has been found for the 
forecasts and intermediate parameter. This is an indication of the correctness of the Vis 
developed , 
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Series A 










Fig. 3.4 Front Panel for series B 







Data ! [Probability 

Series C ^ f 50% 



Fig. 3.5 Front Panel for series C 













38 


Fig. 3.7 Front Panel for series E 
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Fig. 3.8 Front Panel for series F 
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.10 Hierarchy of subVls 
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Fig. 3.11 


Consolidated list of subVIs 
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Fig. 3.14(a) yi for calculating partial autocorrelation 





























Fig. 3.15 VI for calculating initial estimates of autoregressive parameters 






Fig. 3.16(a) VI for calculating the initial estimates of a general ARMA Model 
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Fig. 3.16(b) VI for calculating the initial estimates of a general ARMA Model 
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Fig. 3.16(c) VI for calculating the initial estimates of a general ARMA Model 
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Fig. 3 . 16 (e) VI for calculating the initial estimates of a general ARMA Model 
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Fig. 3.16(f) VI for calculating the initial estimates of a general ARM A Model 
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Fig. 3.17(a) VI to convert (j) to <p 

















Fig. 3.18(a) VI for calculating forecasts at origin t for lead I 






Fig. 3.18(b) VI for c.'ilciilaliiig forecasts at origin t for lead / 
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Fig. 3.18((1) VI for calculating forecasts at origin t for lead / 










Fig. 3.18(e) VI for calculaling forecasts at origin t for lead / 
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Fig. 3.21 VI for calculaliiig a, and S using fig 3.18 and 3.19 
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Fig. 3.22(1) VI for calculating final parameter estimation 
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Fig. 3.23(a) VI for forecasting the final model 
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Fig. 3.23(b) VI for forecasting the final model 
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Fig. 3.23(d) VI for forecasting the final model 
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Fig. 3.23(c) VI for forecasliiig llic final model 
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Fig. 3.23(g) VI for forecasting the final model 
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Fig. 3.23(h) VI for forecasting the final model 
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Fig. 3.23(i) VI for forecasting the final model 























CHAPTER 4 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

During the course of the present work Virtual Instruments for Time Series Analysis have 
been developed. With reference to Figure 1.2, the Vis developed incorporate - 

• Model Identification 

• Parameter Estimation 

• Forecasting 

The programs have been established and tested for data available in literature. 

The remaining tasks of the proposed work are as follows to make a complete Virtual 
Instrumentation package for on-line time series analysis . 

Diagnostic Checking 

The identified and estimated model is to next subjected to diagnostic checks and tests of 
goodness of fit with the observed data. The question of adequacy of the model is to be 
addressed in this stage. The inadequacies need to be defined in statistical terms and the 
model needs to be refined further to accommodate reasonably large body of data spread 
over a reasonable time span. Standard techniques of diagnostic checking involving 
overfitting, autocorrelation check and cumulative periodogram check are proposed to be 
built into the Virtual Instrumentation. 

Noise Considerations 

The problem of model building through a regressive analysis or transfer functions 
becomes complex in the presence of noise in various input/output parameters. Error 
estimation algorithms can be built into the models to check their robustness and to choose 
the model with the smallest possible mean square error. 

On-line Data Acquisition From The Rotor-Rig 

Certain simple instruments like oscilloscopes, FFT Analyzers and Filters were developed 
during the present work in order to check the on-line functioning of the AD card and its 
compatibility with the Lab VIEW software. Vibration data was taken from the rotor-rig 
and displayed successfully on the above Vis. The Time Series Analysis model could not 
be tested on-line due to time constraints. Once the model is complete, by incorporating 
the Noise Considerations described above, it can be made functional to take the data from 
the rotor-rig. 
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